[WIP] Pressure Based Solver - #2812

Proposed Changes

The current work (part of GSoC) provides a working version of a pressure-based algorithm for the incompressible flow solver as an alternative to the existing Density-based solver. Below, the reader may find the algorithm which has been implemented, as well as the current progress of the code and challenges. All the way at the bottom one can find performance comparisons between the DB and PB solvers for some test cases.

Algorithm

A lot of versions of pressure-based algorithms exist, and many different versions can be implemented. Here, we opt for versions of the original SIMPLE/PISO algorithm, it is briefly defined here for clarity.

First, the momentum equations are solved, starting from the previous time step's velocity$\vec{u}^{(0)}$, pressure $p^{(0)}$, and face velocity $\vec{u}_f^{(0)}$. The resulting momentum is the predicted momentum, here its discretized form is shown, as its coefficients are used in the subsequent equations

$A_p\vec{u}_p^{(1)}+\sum_n A_n \vec{u}^{(1)}_n=-\frac{V}{\rho}\nabla p^{(0)}+S_m$

The subsequent momentum is not necessarily incompressible, the pressure correction equation can be derived by rewriting it as follows, using a term often called H by A

$\vec{u}_p^{(1)}=-\frac{\sum_n A_n \vec{u}_n^{(1)}}{A_p}-\frac{V}{\rho A_p}\nabla p+S_m=\frac{H(\vec{u}^{(1)})}{A}-\frac{V}{\rho A_p}\nabla p^{(0)}+S_m$

$\vec{u}_p^{(2)}=\frac{H(\vec{u}^{(1)})}{A}-\frac{V}{\rho A_p}\nabla p^{(1)}+S_m$

Note how a simplification is used here where the HbyA term is neglected. Subtracting these two equations yields the first pressure correction equation for $p'$ as

$\nabla \cdot \left( \frac{V}{\rho A_p}\nabla p'\right)=\nabla \cdot \vec{u}^{(1)}$

Make note that for the divergence here, we require the face mass fluxes, which are computed using Rhie-Chow interpolation to avoid odd-even decoupling. After the equation is solved, using the pressure correction $p'$, the pressure and momentum are corrected according to

$p^{(1)} = p^{(0)} + p' ,\quad u^{(2)} = u^{(1)}+u'.\quad u'= -\frac{V}{\rho A_p}\nabla p'$

So far, this is equal to a pseudo-transient version of the SIMPLE algorithm. This algorithm however suffers from a very tight stability condition on the time-step size. Therefore, multiple pressure corrections can be applied, which for two corrections is originally called the PISO algorithm.

The second pressure correction does not neglect the HbyA term, which then results in the equation

$\nabla \cdot \left( \frac{V}{\rho A_p}\nabla p'\right)=\nabla \cdot \vec{u}^{(2)}+\nabla\cdot\left(\frac{H(\vec{u}')}{A}\right)$

And the new correction equations are defined as

$p^{(2)} = p^{(1)} + p' ,\quad u^{(3)} = u^{(2)}+u'.\quad u'= \frac{H(\vec{u}')}{A}-\frac{V}{\rho A_p}\nabla p'$

Note that HbyA here uses the previous velocity correction and is thus the same quantity as the one used in the second pressure correction equation. Later pressure correction equations follow analogously.

Progress:

• Pressure-based solver added as alternative to density-based solver for the incompressible flow equations.
• The pressure-based solver has only been tested for constant density cases for basic Navier-Stokes and Euler flow.
• The pressure-based solver is implemented based on a pseudo time-stepping approach to remain consistent with the other solvers in SU2. The pressure-based solver is currently set to the SIMPLE algorithm by default, with options for SIMPLEC and PISO available.
• The Poisson solver is a major bottleneck in the computation speed. To account for this, an option was added to use a different linear solver and preconditioner for the poisson solver.
• For details on the implementation of the algorithm and the responsibility distribution please see the file CPBFluidIteration.cpp.
• A tutorial showcasing the different options of the PB solver for the lid driven cavity flow problem has been added, please see Lid driven cavity flow tutorial Tutorials#86 and Add Lid Driven Cavity Flow tutorial for pressure-based solver su2code.github.io#218

Issues

Performance:

• The Poisson solver can sometimes struggle a lot due to high Reynolds numbers and fine meshes, and thus require a ridiculous number of iterations to converge reasonably. Possible fixes include adding multigrid support or a DIC preconditioner (far less efficient). Multigrid support is tricky as SU2 currently only considers multigrid for the main (flow) solver and not for auxiliary solvers.

• Convergence issues with RANS (SA and SST) on fine meshes with high Reynolds numbers. Tests have shown that cases such as flow over a flat plate converges fine. However, external aerodynamic cases such as the naca0012 RANS test case do not converge well at all. The convergence does slightly improve when we switch out the mesh for a more uniform unstructured mesh without large aspect ratio cells in the wake of the airfoil, although this only slightly helps. The flat plate turbulence test case also uses large aspect ratio cells so this is not the sole issue. The RANS solver also often requires many iterations of the Poisson solver to converge reasonably, this is however not the reason for the lack of convergence.

• Periodic boundary conditions have not been implemented/tested at all as of yet.

• Any code related to adjoints has not been considered at all either.

Code:

• Parallelization with OMP gives wrong results, MPI however does work as expected.

TODO list

• Fix issues mentioned above (future work)
• Fix code so it can be merged

Related Work

This work is based on earlier attempts by Nitish Anand (2024) and Akshay Koodly (2021), see feature branches feature_PBFlow_V8 and feature_Pressure_based respectively. Also see PR #2210

PR Checklist

• I am submitting my contribution to the develop branch.
• My contribution generates no new compiler warnings (try with --warnlevel=3 when using meson).
• My contribution is commented and consistent with SU2 style (https://su2code.github.io/docs_v7/Style-Guide/).
• I used the pre-commit hook to prevent dirty commits and used pre-commit run --all to format old commits.
• I have added a test case that demonstrates my contribution, if necessary.
• I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp), if necessary.

Result showcase

Inviscid Hydrofoil

Convergence history of the inviscid flow around a hydrofoil at a 5 degree aoa.

The pressure coefficient along the surface of the hydrofoil at a 5 degree aoa and the corresponding lift coefficients, X-FOIL predicts C_L=0.6.

Lid Driven Cavity

Convergence history of the lid driven cavity problem, note that CFL=60 is the highest stable CFL for the PB solver, whereas the DB solver does not have this CFL related stability issue.

Flatplate RANS

The skin friction coefficient for turbulent flow over a (rough) flat plate with SA.